%%% proj_meta_testing

enu_ind=[1 2 4 5 11 12 15 16];
ctr_ind=setdiff([1:16],enu_ind);

t_points=[-48 -24 12 24 48 72];

xc_enu=[];
tl_enu=[];
cnt=0;
for siteID=enu_ind
    
    cnt=cnt+1;
    for tpoint=1:6
        grat_onset_frame=find(round(proj_meta(siteID).rd(1,tpoint).stim_id)>0,1,'first');
        frames_to_use=round(proj_meta(siteID).rd(1,tpoint).ps_id/10)==1;
        frames_to_use(grat_onset_frame:end)=0;
        
        all_act=[proj_meta(siteID).rd(1,tpoint).act;proj_meta(siteID).rd(2,tpoint).act;proj_meta(siteID).rd(3,tpoint).act;proj_meta(siteID).rd(4,tpoint).act];
        mean_act=mean(all_act);
        tmp_a=mean_act(frames_to_use);
        tmp_b=proj_meta(siteID).rd(1,tpoint).velM_smoothed(frames_to_use);
        tmp_a=tmp_a-mean(tmp_a);
        tmp_b=tmp_b-mean(tmp_b);
        tmp_a=tmp_a/max(tmp_a);
        tmp_b=tmp_b/max(tmp_b);
%         figure;plot(tmp_a);hold on;plot(tmp_b,'r')
        tmp=xcorr(tmp_a,tmp_b,'coeff');
        [tmp_val,tmp_pos]=max(tmp);
        xc_enu(cnt,tpoint)=tmp_val;
        tl_enu(cnt,tpoint)=tmp_pos;
%             figure;
%             imagesc(all_act)
%             set(gca,'clim',[1 4])
%             figure; hold on
%             plot(mean(all_act))
%             plot(proj_meta(siteID).rd(1,tpoint).velM_smoothed*10+1,'r')
%             plot(proj_meta(siteID).rd(1,tpoint).ps_id/20+0.9,'g')
        
        
    end
end

xc_ctr=[];
tl_ctr=[];
cnt=0;
for siteID=ctr_ind
    cnt=cnt+1;
    for tpoint=1:6
        grat_onset_frame=find(round(proj_meta(siteID).rd(1,tpoint).stim_id)>0,1,'first');
        frames_to_use=round(proj_meta(siteID).rd(1,tpoint).ps_id/10)==1;
        frames_to_use(grat_onset_frame:end)=0;
        
        all_act=[proj_meta(siteID).rd(1,tpoint).act;proj_meta(siteID).rd(2,tpoint).act;proj_meta(siteID).rd(3,tpoint).act;proj_meta(siteID).rd(4,tpoint).act];
        mean_act=mean(all_act);
        tmp_a=mean_act(frames_to_use);
        tmp_b=proj_meta(siteID).rd(1,tpoint).velM_smoothed(frames_to_use);
        tmp_a=tmp_a-mean(tmp_a);
        tmp_b=tmp_b-mean(tmp_b);
        tmp_a=tmp_a/max(tmp_a);
        tmp_b=tmp_b/max(tmp_b);
        tmp=xcorr(tmp_a,tmp_b,'coeff');
        [tmp_val,tmp_pos]=max(tmp);
        xc_ctr(cnt,tpoint)=tmp_val;
        tl_ctr(cnt,tpoint)=tmp_pos;
    end
end

figure; hold on
plot(t_points(1:2),mean(xc_ctr(:,1:2)),':k','linewidth',2)
plot(t_points(1:2),mean(xc_enu(:,1:2)),':r','linewidth',2)
plot(t_points(3:6),mean(xc_ctr(:,3:6)),':k','linewidth',2)
plot(t_points(3:6),mean(xc_enu(:,3:6)),':r','linewidth',2)

for ind=1:length(t_points)
    plot([t_points(ind) t_points(ind)],mean(xc_ctr(:,ind))+[-1 1]*std(xc_ctr(:,ind))/sqrt(size(xc_ctr,1)),'k');
    plot([t_points(ind) t_points(ind)],mean(xc_enu(:,ind))+[-1 1]*std(xc_enu(:,ind))/sqrt(size(xc_enu,1)),'r');
end

plot(t_points,mean(xc_ctr),'ko','markersize',8,'linewidth',2)
plot(t_points,mean(xc_enu),'ro','markersize',8,'linewidth',2)

set(gca,'xtick',[t_points(1:2) 0 t_points(3:6)])
title('correlation of neural activity with locomotion')
legend('ctr','enu')



avg_act=[];
cnt=0;
for siteID=enu_ind
    for pl=1:4
        for tpoint=1:6
            grat_onset_frame=find(round(proj_meta(siteID).rd(1,tpoint).stim_id)>0,1,'first');
            frames_light=round(proj_meta(siteID).rd(1,tpoint).ps_id/10)==1;
            frames_dark=round(proj_meta(siteID).rd(1,tpoint).ps_id/10)==0;
            frames_light(grat_onset_frame:end)=0;
            frames_dark(grat_onset_frame:end)=0;
            frames_running=proj_meta(siteID).rd(pl,tpoint).velP_smoothed>0.005;
            
            
            avg_act(tpoint,cnt+1:cnt+size(proj_meta(siteID).rd(pl,tpoint).act,1),1)=mean(proj_meta(siteID).rd(pl,tpoint).act(:,frames_light&~frames_running)');
            avg_act(tpoint,cnt+1:cnt+size(proj_meta(siteID).rd(pl,tpoint).act,1),2)=mean(proj_meta(siteID).rd(pl,tpoint).act(:,frames_light&frames_running)');
            avg_act(tpoint,cnt+1:cnt+size(proj_meta(siteID).rd(pl,tpoint).act,1),3)=mean(proj_meta(siteID).rd(pl,tpoint).act(:,frames_dark&~frames_running)');
            avg_act(tpoint,cnt+1:cnt+size(proj_meta(siteID).rd(pl,tpoint).act,1),4)=mean(proj_meta(siteID).rd(pl,tpoint).act(:,frames_dark&frames_running)');
            
        end
        
        cnt=cnt+size(proj_meta(siteID).rd(pl,tpoint).act,1);
    end
end

act_cells=(squeeze(max(max(avg_act,[],1),[],3))>1.5);
avg_aa=avg_act(:,act_cells,:);
for sort_cond=1:4
    [~,sid]=sort((mean(avg_aa(1:2,:,sort_cond))-mean(avg_aa(3:6,:,sort_cond)))./(mean(avg_aa(1:2,:,sort_cond))+mean(avg_aa(3:6,:,sort_cond))));
    figure;
    for knd=1:4
        subplot(2,2,knd)
        imagesc(avg_aa(:,sid,knd)')
        set(gca,'clim',[1 1.7])
    end
end

clrs='rgbc';
figure;hold on
for ind=1:4
    plot(mean(avg_aa(:,:,ind)'),'color',clrs(ind));
end
title('Avg ENU activity')
legend('sitting light','running light','sitting dark','running dark')

tmp_pre=squeeze(mean(avg_aa(1:2,:,:),1));
figure; hold on
plot(tmp_pre(:,3),tmp_pre(:,4),'.')
plot(tmp_pre(:,2),tmp_pre(:,4),'g.')
plot(tmp_pre(:,2),tmp_pre(:,1),'r.')

tmp_post=squeeze(mean(avg_aa(3:6,:,:),1));
figure; hold on
plot(tmp_post(:,3),tmp_post(:,4),'.')
plot(tmp_post(:,2),tmp_post(:,4),'g.')
plot(tmp_post(:,2),tmp_post(:,1),'r.')


figure;
plot(mean(tmp_pre(:,1:2),2)-mean(tmp_pre(:,3:4),2),mean(tmp_post(:,1:2),2)-mean(tmp_post(:,3:4),2),'g.')


fp=[];
for ind=1:4
    fp(ind,:)=mean(avg_aa(1:2,:,ind),1);
end


num_clust=4;
idx=kmeans(fp',num_clust);

for ind=1:num_clust
    figure;
    for knd=1:4
        subplot(2,2,knd)
        imagesc(avg_aa(:,idx==ind,knd)')
        set(gca,'clim',[1 1.7])
    end
    figure;
    imagesc(fp(:,idx==ind))
    set(gca,'clim',[1 1.7])
end

for ind=1:num_clust
    figure;plot(squeeze(median(avg_aa(:,idx==ind,:),2)))
end














for site_id=1:size(proj_meta,2)
    nbr_time_points_site(site_id)=size(proj_meta(site_id).rd,2);
end

ccfs=nan(max(nbr_time_points_site),max(nbr_time_points_site),size(proj_meta,2));
ccfs_l=ccfs; ccfs_d=ccfs; ccfs_g=ccfs;



fb_da_cc=nan(max(nbr_time_points_site),size(proj_meta,2));
for site_id=1:size(proj_meta,2)
    act_cur_site=[];
    light_act=[]; dark_act=[]; grat_act=[];
    for time_point=1:size(proj_meta(site_id).rd,2)
        cur_act=[proj_meta(site_id).rd(1,time_point).act;proj_meta(site_id).rd(2,time_point).act;proj_meta(site_id).rd(3,time_point).act;proj_meta(site_id).rd(4,time_point).act];
        
        act_cur_site(:,time_point)=mean(cur_act,2);
        
        gratings_start=find(abs(round(proj_meta(site_id).rd(1,time_point).stim_id))>0,1,'first');
        
        dark_frames=round(proj_meta(site_id).rd(1,time_point).ps_id)==0;
        light_frames=round(proj_meta(site_id).rd(1,time_point).ps_id)>0;
        light_frames(gratings_start:end)=0;
        grat_frames=zeros(length(light_frames),1);
        grat_frames(gratings_start:end)=1;
        grat_frames=logical(grat_frames);
        
        light_act(:,time_point)=mean(cur_act(:,light_frames),2);
        dark_act(:,time_point)=mean(cur_act(:,dark_frames),2);
        grat_act(:,time_point)=mean(cur_act(:,grat_frames),2);
    end
    for time_point=1:size(proj_meta(site_id).rd,2)
        tmp=corrcoef(light_act(:,time_point),dark_act(:,time_point));
        ld_cc(time_point,site_id)=tmp(2);
        tmp=corrcoef(light_act(:,time_point),grat_act(:,time_point));
        lg_cc(time_point,site_id)=tmp(2);
        tmp=corrcoef(grat_act(:,time_point),dark_act(:,time_point));
        gd_cc(time_point,site_id)=tmp(2);
        
        
        tmp=corrcoef(light_act(:,time_point),mean(dark_act(:,1:2),2));
        ld1_cc(time_point,site_id)=tmp(2);
        tmp=corrcoef(light_act(:,time_point),mean(light_act(:,1:2),2));
        ll1_cc(time_point,site_id)=tmp(2);
        
        tmp=corrcoef(dark_act(:,time_point),mean(dark_act(:,1:2),2));
        dd1_cc(time_point,site_id)=tmp(2);
        tmp=corrcoef(dark_act(:,time_point),mean(light_act(:,1:2),2));
        dl1_cc(time_point,site_id)=tmp(2);
    end

    tmp=corrcoef(act_cur_site);
    ccfs(1:size(tmp,1),1:size(tmp,2),site_id)=tmp;
    
    tmp=corrcoef(light_act);
    ccfs_l(1:size(tmp,1),1:size(tmp,2),site_id)=tmp;
    
    tmp=corrcoef(dark_act);
    ccfs_d(1:size(tmp,1),1:size(tmp,2),site_id)=tmp;
    
    tmp=corrcoef(grat_act);
    ccfs_g(1:size(tmp,1),1:size(tmp,2),site_id)=tmp;

end

figure;imagesc(mean(ccfs(:,:,enu_ind),3))
set(gca,'clim',[0 1])
figure;imagesc(mean(ccfs(:,:,ctr_ind),3))
set(gca,'clim',[0 1])


figure;imagesc(mean(ccfs_l(:,:,enu_ind),3))
set(gca,'clim',[0 1])
title('ENU corr light activity')
figure;imagesc(mean(ccfs_l(:,:,ctr_ind),3))
set(gca,'clim',[0 1])
title('CTR corr light activity')
figure;imagesc(mean(ccfs_d(:,:,enu_ind),3))
set(gca,'clim',[0 1])
title('ENU corr dark activity')
figure;imagesc(mean(ccfs_d(:,:,ctr_ind),3))
set(gca,'clim',[0 1])
title('CTR corr dark activity')
figure;imagesc(mean(ccfs_g(:,:,enu_ind),3))
set(gca,'clim',[0 1])
title('ENU corr grating activity')
figure;imagesc(mean(ccfs_g(:,:,ctr_ind),3))
set(gca,'clim',[0 1])
title('CTR corr grating activity')


figure; hold on
plot(ld_cc(:,ctr_ind))
plot(mean(ld_cc(:,ctr_ind)'),'k','linewidth',2)


figure; hold on
plot(ld_cc(:,enu_ind))
plot(mean(ld_cc(:,enu_ind)'),'k','linewidth',2)

figure;hold on
plot(t_points(1:2),mean(ld_cc(1:2,enu_ind)'),'r','linewidth',2)
plot(t_points(1:2),mean(ld_cc(1:2,ctr_ind)'),'k','linewidth',2)
plot(t_points(3:6),mean(ld_cc(3:6,enu_ind)'),'r','linewidth',2)
plot(t_points(3:6),mean(ld_cc(3:6,ctr_ind)'),'k','linewidth',2)
plot(t_points,mean(ld_cc(:,enu_ind)'),'ro','linewidth',2)
plot(t_points,mean(ld_cc(:,ctr_ind)'),'ko','linewidth',2)
title('correlation of activity patterns between light and dark running')



figure;hold on;
plot(mean(ld1_cc(:,enu_ind),2),'r')
plot(mean(ll1_cc(:,enu_ind),2),'r')
plot(mean(ld1_cc(:,ctr_ind),2),'k')
plot(mean(ll1_cc(:,ctr_ind),2),'k')

figure;hold on;
plot(mean(dd1_cc(:,enu_ind),2),'r')
plot(mean(dl1_cc(:,enu_ind),2),'r')
plot(mean(dd1_cc(:,ctr_ind),2),'k')
plot(mean(dl1_cc(:,ctr_ind),2),'k')





%%%%%%%%%%%%%%%%%%%%%%% PCA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

all_act=[];

site_id=12;

for tp=1:6
    y_start=size(all_act,2);
    x_start=0;
    l_ind=(proj_meta(site_id).rd(zl,tp).ps_id<1);
    l_ind(grat_onset:end)=logical(0);
    
    for zl=1:4
        tmp=proj_meta(site_id).rd(zl,tp).act;
        grat_onset=find(proj_meta(site_id).rd(zl,tp).stim_id>0.5,1,'first');
        all_act(x_start+1:x_start+size(tmp,1),y_start+1:y_start+sum(l_ind))=tmp(:,l_ind);
        x_start=x_start+size(tmp,1);
    end
end

pc=princomp(all_act');
figure;imagesc(all_act)
figure;imagesc(pc)

pc10=pc(:,1:10);

tmp=all_act'*pc10;

figure;imagesc(tmp')
























